* Analysis on ADM1 Level

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
********************
* Data Preparation *
********************
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
set maxvar 32767
**# Bookmark #2
use "$data\Merge\IDA_Aid_GED_ADM1.dta", clear

* keep only African dataset
keep if unregion2=="Africa"
// drop small economies and islands
drop if d_smallcountry ==1
* Generate continental identifiers
egen continent_num=group(unregion2)


* Generate Trend Variables:
sort ID_adm1_num transaction_year
bysort ID_adm1_num: egen mean_temp=mean(temp)
egen trend=group(transaction_year)
gen trend_sq=trend*trend
gen trend_cub=trend*trend*trend

* Aid data
gen lnaid=ln(WBAID_ADM1_LOC+0.01)
label var lnaid  "Log of location weighted World Bank Aid per capita +0.01"
gen lnaid_c=ln(CAID_ADM1_LOC+0.01)
label var lnaid_c  "Log of location weighted Chinese Aid per capita +0.01"
gen lnaid_pop=ln(WBAID_ADM1_POP+0.01)
label var lnaid_pop  "Log of population weighted World Bank Aid per capita +0.01"
gen lnaid_pop_c=ln(CAID_ADM1_POP+0.01)
label var lnaid_pop_c  "Log of population weighted Chinese Aid per capita +0.01"

* Conflict Data
* Rename SCAD variables for consistency with other labeling
renvars D_pro_gov_violence D_pro_gov_vio_leth D_pro_gov_vio_nleth D_demo D_riot D_strike D_total_violence D_total_rds / D_pro_gov_v_adm1 D_pro_gov_v_l_adm1 D_pro_gov_v_nl_adm1 D_demo_adm1 D_riot_adm1 D_strike_adm1 D_tot_v_adm1 D_tot_rds_adm1
* Rescale Dummies for easier interpretation
foreach var in D_pro_gov_v_adm1 D_pro_gov_v_l_adm1 D_pro_gov_v_nl_adm1 D_demo_adm1 D_riot_adm1 D_strike_adm1 D_tot_v_adm1 D_tot_rds_adm1 D_repress_all D_repress_nl D_repress_l {
	replace `var'=`var'*100
}
gen lnbrd_adm1=ln(brd_adm1+0.01)
* Genreate new intensity thresholds
gen inten1_adm1=0
replace inten1_adm1=100 if brd_adm1>=5
gen inten2_adm1=0
replace inten2_adm1=100 if brd_adm1>=25

* Prepare Ending and Onset Indicators
sort ID_adm1_num transaction_year
xtset ID_adm1_num transaction_year
gen d1inten1_adm1=d1.inten1_adm1
* Onset and Ending on regional level
gen onset_adm1=0
replace onset_adm1=100 if d1inten1_adm1==100
gen ending_adm1=0
replace ending_adm1=100 if d1inten1_adm1==-100
* Onset and Ending on country level
bysort state transaction_year: egen conflict_country=max(inten1_adm1)
sort ID_adm1_num transaction_year
gen d1conflict_country=d1.conflict_country
gen onset_country=0
replace onset_country=100 if d1conflict_country==100
gen ending_country=0
replace ending_country=100 if d1conflict_country==-100

* Generate Latitude and Longitude for time-interaction
gen lon=abs(xcoord)
gen lat=abs(ycoord)

*** Prepare Rolling IV - WB Aid***
*2 Years
xtset ID_adm1_num transaction_year
gen received=0
replace received=1 if WBAID_ADM1_LOC!=0 & WBAID_ADM1_LOC!=. 
gen eligible=1
sort ID_adm1_num transaction_year
bysort ID_adm1_num: gen eligible_sum=sum(eligible)
bysort ID_adm1_num: gen received_sum=sum(received)
bysort ID_adm1_num: egen eligible_tot=total(eligible)
bysort ID_adm1_num: egen received_tot=total(received)
gen prob_cum=received_sum/eligible_sum
gen prob_cons=received_tot/eligible_tot
gen prob_rolling2=(l1.received+l0.received)/2 if transaction_year>=1996
gen prob_rolling3=(l2.received+l1.received+l0.received)/3 if transaction_year>=1997
gen prob_rolling4=(l3.received+l2.received+l1.received+l0.received)/4 if transaction_year>=1998
drop eligible* received*


*** Prepare Rolling IV - Chinese Aid***
xtset ID_adm1_num transaction_year
*2 Years
gen received=0
replace received=1 if CAID_ADM1_LOC!=0 & CAID_ADM1_LOC!=. & transaction_year>=2000 & transaction_year<=2012
gen eligible=1 if transaction_year>=2000 & transaction_year<=2012
sort ID_adm1_num transaction_year
bysort ID_adm1_num: gen eligible_sum=sum(eligible) if transaction_year>=2000 & transaction_year<=2012
bysort ID_adm1_num: gen received_sum=sum(received) if transaction_year>=2000 & transaction_year<=2012
bysort ID_adm1_num: egen eligible_tot=total(eligible) if transaction_year>=2000 & transaction_year<=2012
bysort ID_adm1_num: egen received_tot=total(received) if transaction_year>=2000 & transaction_year<=2012
gen prob_cum_c=received_sum/eligible_sum if transaction_year>=2000 & transaction_year<=2012
gen prob_cons_c=received_tot/eligible_tot if transaction_year>=2000 & transaction_year<=2012
gen prob_rolling_c2=(l1.received+l0.received)/2 if transaction_year>=2001
gen prob_rolling_c3=(l2.received+l1.received+l0.received)/3 if transaction_year>=2002
gen prob_rolling_c4=(l3.received+l2.received+l1.received+l0.received)/4 if transaction_year>=2003
drop eligible* received*

* Detrend Chinese steel production:
* Lennart 23.03.2018: Why should we add 0.01 as steel production is always positive & larger zero between 2000-2012
*create ln of steel variable
gen steel_prod_ln=ln(steel_prod+0.01)
capture gen constant=1
reg steel_prod_ln trend constant, nocons
gen steel_prod_ln_detrend = steel_prod_ln-_b[trend]*trend-_b[constant]

*Detrend detrended Factor loading of other Chinese raw materials (10.09.2021: Previously erroneous code)
*capture gen constant=1
*reg factor1 trend constant, nocons
*gen factor1_detrend = factor1-_b[trend]*trend-_b[constant]

* Prepare IBRD variable
gen temp_IBRD = IBRD_ADM1_LOC
replace temp_IBRD=0 if IBRD_ADM1_LOC<=0
gen ln_IBRD_ADM1_LOC=ln(temp_IBRD+0.01)
drop temp_IBRD

* Create IDA_position fluctuation cleaned from general trends
reg IDA_position trend
predict IDA_position_cleaned, resid

* Create running average of IDA & IBRD position (e.g., for aid in t IDA_position in t determines the first half of the transaction_year and aid in t+1 determines the second half of the transaction_year)
sort ID_adm1_num transaction_year
gen IDA_position_run=(IDA_position+f1.IDA_position)/2
gen IBRD_EtL_Ratio_run=(IBRD_EtL_Ratio+f1.IBRD_EtL_Ratio)/2
* For out of sample values refer to IDA_IBRD values from Brad Parks (see Data/Aid/IDA_IBRD_Liquidity.xlsx)
*replace IBRD_EtL_Ratio_run=(0.268+0.2698)/2 if transaction_year==2012
*replace IDA_position_run=(0.788+0.807)/2 if transaction_year==2012
replace IBRD_EtL_Ratio_run=(0.217+0.2142)/2 if transaction_year==1995
replace IDA_position_run=(1.106+0.981)/2 if transaction_year==1995

* Detrend IDA Position
capture gen constant=1
reg IDA_position_run trend constant, nocons
gen IDA_position_run_detrend = IDA_position_run-_b[trend]*trend-_b[constant]

* Set Chinese aid to missing for non-African states and year<2000 & years>2012
replace lnaid_c=. if unregion2!="Africa" | transaction_year<2000 | transaction_year>2012

* Prepare first difference analysis
gen d1lnaid=d1.lnaid
gen d1lnaid_c=d1.lnaid_c
gen d1IDA_position_run=d1.IDA_position_run
gen d1steel_prod_ln=d1.steel_prod_ln
gen d1factor1=d1.factor1

* Generate identifiers for clusters
egen fcy=group(state transaction_year)
egen fc=group(state)
egen fr=group(ID_adm1_num)

* Labeling
label var lnbrd_adm1  "log of Battle-related deaths + 0.01"
label var ln_IBRD_ADM1_LOC "log of location weighted IBRD funding"
label var inten1_adm1 "dichotomous indicator (0/100) for for 5>=battle deaths"
label var inten2_adm1 "dichotomous indicator (0/100) for for 25>=battle deaths "
label var d1inten1_adm1 "first difference of inten1_adm1"
label var onset_adm1 "Dummy for an intensity 1 conflict starting at the regional level with no prior conflict"
label var ending_adm1 "Dummy for an intensity 1 conflict ending at the regional level with no prior conflict"
label var conflict_country "dichotomous indicator (0/100) for at least 5>=battle deaths in any ADM region"
label var d1conflict_country "first difference of conflict_country"
label var onset_country "Dummy for intensity 1 conflict starting at the regional level for country with no prior conflict"
label var ending_country "Dummy for all regional intensity 1 conflict ending in the country"
label var lon "Longitude"
label var lat "Latitude"
label var prob_cum "IDA: Cumulative/time variant probability"
label var prob_rolling2 "IDA: Cumulative/time variant probability for current + 1 previous year"
label var prob_rolling3 "IDA: Cumulative/time variant probability for current + 2 previous years"
label var prob_rolling4 "IDA: Cumulative/time variant probability for current + 3 previous years"
label var prob_cons "IDA: Constant/non-time variant probability"
label var prob_cum_c "China: Cumulative/time variant probability"
label var prob_cons_c "China: Constant/non-time variant probability"
label var prob_rolling_c2 "China: Cumulative/time variant probability for current + 1 previous year"
label var prob_rolling_c3 "China: Cumulative/time variant probability for current + 2 previous years"
label var prob_rolling_c4 "China: Cumulative/time variant probability for current + 3 previous years"
label var steel_prod_ln "Log of Chinese steel production"
label var steel_prod_ln_detrend "Log of Chinese steel production without linear trend"
label var IDA_position_run "Two year moving average of IDA position to account for fiscal year"
label var IBRD_EtL_Ratio_run "Two year moving average of IBRD funding ratio to account for fiscal year"
label var IDA_position_run_detrend "IDA_position_run without linear trend"
label var d1IDA_position_run "first difference of IDA_position_run"
label var d1steel_prod_ln "first difference of steel_prod_ln"
label var d1factor1 "first difference of factor1"
label var d1lnaid "first difference of lnaid"
label var d1lnaid_c  "first difference of lnaid_c"
label var trend_sq "Quadratic trend"
label var trend_cub "Cubic trend"
label var mean_temp "Constant/non-time variant regional temperature"
label var factor1 "Scores for factor 1 (detrended commodity factor inputs)"
label var factor1_levs "Scores for factor 1 levels (non-detrended commodity factor inputs)"

***************
* Globals for *
* Regressions *
***************
xtset ID_adm1_num transaction_year
* All controls
global control1 droughtend_spi droughtstart_spi temp prec_gpcc 

global control2 droughtend_spi droughtstart_spi temp prec_gpcc c.ttime_mean#i.transaction_year  c.landarea#i.transaction_year  c.dist_capital#i.transaction_year c.mean_temp#i.transaction_year c.elevation_std#i.transaction_year c.Dborder#i.transaction_year c.Driver#i.transaction_year c.Dlake#i.transaction_year c.Docean#i.transaction_year c.lat#i.transaction_year c.lon#i.transaction_year
*Labeling
global ln "$\ln" 
xtset ID_adm1_num transaction_year

******************Standardize and transfrom variables ROSSELLA
cap drop std_lnaid std_lnaid_c 
egen std_lnaid=std(lnaid), mean(0)  std(1)
egen std_lnaid_c=std(lnaid_c), mean(0)  std(1)
gen asinhimp_glob=asinh(Imports_global)
gen asinhfdi_glob=asinh(FDI_global)
gen asinhbrd_global=asinh(brd_global)
gen asinhgdp_glob=asinh(GDP_global)
bysort ADM0 transaction_year: egen brd_adm0=total(brd_adm1)
gen asinhbrd_adm0=asinh(brd_adm1)
*************************************
*Variables for future regressions
gen nat_res_s=.
replace nat_res_s=0 if !missing(petroleum_s, diamsec_s, diamprim_s, gem_s, goldvein_s, goldsurface_s, goldplacer_s, gem_s)
replace nat_res_s=1 if petroleum_s==1 | diamsec_s==1 | diamprim_s==1 | gem_s==1 | goldvein_s==1 | goldsurface_s==1 | goldplacer_s==1 | gem_s==1 
foreach var in chinese_fdi chinese_imports global_fdi global_imports {
gen asinh_`var'=asinh(`var')
}
* Generate new stable probability from third year of panel
gen aux_wb=.
replace aux_wb=prob_cum if transaction_year==1998
bysort ID_adm1: egen prob_cum_98=max(aux_wb)
gen aux_c=.
replace aux_c=prob_cum_c if transaction_year==2003
bysort ID_adm1: egen prob_cum_c_03=max(aux_c)
drop aux*



gen aux_wb=.
replace aux_wb=prob_cum if transaction_year==2003
bysort ID_adm1: egen prob_cum_03=max(aux_wb)
gen aux_c=.
replace aux_c=prob_cum_c if transaction_year==2006
bysort ID_adm1: egen prob_cum_c_06=max(aux_c)
drop aux*

**# Bookmark #17
*Sort
sort ID_adm1_num transaction_year
xtset ID_adm1_num transaction_year

save "$data\working_file.dta",replace
**# Bookmark #1


***********************
* IV Robustness Tests *
***********************
use "$data\working_file.dta",clear
**********
* Graphs *
**********

* Label IDA_position
label var IDA_position  "IDA Position"
label var IDA_position_run "Running average of IDA Position"


* Generate graphs of (un)common trends

bysort transaction_year: egen mean_brd_low=mean(brd_adm1) if prob_cons<0.5
bysort transaction_year: egen mean_brd_high=mean(brd_adm1) if prob_cons>0.5
twoway line  IDA_position_run transaction_year, lpattern(solid) lwidth(thick) yaxis(2) yscale(axis(2) range(0.6(0.1)1.2)) yla(0.6(0.1)1.2) || line mean_brd_low transaction_year, yaxis(1) yscale(alt) yscale(alt axis(2)) lcolor(blue)  lpattern(shortdash) lwidth(thick) || line mean_brd_high transaction_year,  lpattern(longdash) lwidth(thick)  lcolor(orange) yaxis(1) yscale(alt) yscale(alt axis(2))  xtitle("Year") ytitle("BRD per Region") name(gr5, replace)   bgcolor(white)  graphregion( color(white)) legend(label(1 "IDA Position") label(2 "BRD Low Prob.") label(3 "BRD High Prob.") cols(1) size(3))
drop mean_brd_low mean_brd_high 
graph export "$graphs\low_vs_high_prob_cons_5_WB.png", as(png) replace
graph close



*****************
*****************
*****************
* Bluhm et al. Factor 1
xtset ID_adm1_num transaction_year

bysort transaction_year: egen mean_inten1_low=mean(inten1_adm1) if prob_cons_c<0.05 & transaction_year>=2000 & transaction_year<=2012
bysort transaction_year: egen mean_inten1_high=mean(inten1_adm1) if prob_cons_c>0.05 & transaction_year>=2000 & transaction_year<=2012
bysort transaction_year: egen mean_lnaid_c_low=mean(lnaid_c) if prob_cons_c<0.05 & transaction_year>=2000 & transaction_year<=2012
bysort transaction_year: egen mean_lnaid_c_high=mean(lnaid_c) if prob_cons_c>0.05 & transaction_year>=2000 & transaction_year<=2012

twoway  line factor1 transaction_year if transaction_year>=2000 & transaction_year<=2012, yaxis(2) lpattern(solid) lwidth(thick)  yscale( axis(2) range(-1.8(0.2)2)) yla(-1.8(0.2)2) || line mean_inten1_low transaction_year if transaction_year>=2000 & transaction_year<=2012, yaxis(1)  lcolor(blue) lwidth(thick) yscale(alt axis(1) range(0(2)22)) lpattern(longdash) lwidth(thick)    || line mean_inten1_high transaction_year if transaction_year>=2000 & transaction_year<=2012,  lpattern(shortdash) lcolor(orange) lwidth(thick) yaxis(1) yscale(alt axis(1) range(0(2)22)) yla(0(2)22)   xtitle("Year", margin(medsmall)) ytitle("Conflict Prob. per Region", axis(1) margin(medsmall)) ytitle("Chinese Commodities Production", axis(2) margin(medsmall))   legend(label(1 "Chinese Commodities") label(2 "Conflict Prob. Low Aid Prob.") label(3 "Conflict Prob. High Aid Prob.")  size(3) cols(1)) name(gr05, replace)   graphregion( color(white) ) bgcolor(white)
graph export "$graphs\low_vs_high_prob_cons_c_05_China1_conflict_nodetrend.png", as(png) replace
graph close
sum mean_inten1_low mean_inten1_high
drop mean_inten1_low mean_inten1_high  mean_lnaid_c_high mean_lnaid_c_low



* Fabricated trends
* IDA
gen mean_brd_high_fabric_wb=1.1-(transaction_year-1995)*0.02
gen mean_brd_low_fabric_wb=0.7
twoway  line   IDA_position_run transaction_year, yaxis(1) yscale(axis(1) range(0.7(0.1)1.1)) yla(0.7(0.1)1.1) lpattern(solid) lwidth(thick)  || line mean_brd_low_fabric_wb transaction_year, yaxis(2)  yscale(alt axis(2) range(0.65(0.1)1.2)) yla(0.65(0.1)1.2) lcolor(blue) lpattern(longdash) lwidth(thick)  || line mean_brd_high_fabric_wb transaction_year,   lpattern(shortdash) lwidth(thick) lcolor(orange) yaxis(2)  yscale(alt axis(2) range(0.65(0.1)1.2)) yla(0.65(0.1)1.2)  xtitle("Year", margin(medsmall)) ytitle("Conflict Prob. per Region", margin(medsmall) axis(1)) ytitle("IDA Position", margin(medsmall) axis(2)) name(gr`t', replace)   legend(label(1 "IDA Position") label(2 "Conflict Prob. Low Aid Prob. (fabric.)") label(3 "Conflict Prob. Low Aid Prob. (fabric.)") cols(1) size(3)) graphregion( color(white) ) bgcolor(white)
graph export "$graphs\low_vs_high_brd_fabricated_WB_conflict.png", as(png) replace



gen mean_brd_high_fabric_china1=-0.08*(((transaction_year-2000)-7)^2)+3.2 if transaction_year>=2000
gen mean_brd_low_fabric_china1= 2 if transaction_year>=2000
* Bluhm et al. (2019) Factor 1
twoway  line factor1 transaction_year if transaction_year>=2000 & transaction_year<=2012, yaxis(2) lpattern(solid) lwidth(thick)  yscale( axis(2) range(-0.2 0.225)) yla(-0.2(0.1)0.225) || line mean_brd_low_fabric_china1 transaction_year if transaction_year>=2000 & transaction_year<=2012,   yaxis(1) yscale(alt) lcolor(blue) lwidth(thick) yscale(alt axis(1) range(-2 3)) lpattern(longdash) ///
 || line mean_brd_high_fabric_china1 transaction_year if transaction_year>=2000 & transaction_year<=2012,  lpattern(shortdash)  lwidth(thick) ///
  yaxis(1) yscale(alt) yscale(alt axis(1) range(-2 3)) yla(-2 (1) 3)  xtitle("Year", margin(medsmall)) ytitle("Conflict Probability per Region", axis(1) margin(medsmall)) ytitle("Chinese Commodities Production", axis(2) margin(medsmall))  lcolor(orange)  legend(label(1 "Chinese Commodities") label(2 "Conflict Prob. Low Aid Prob. (fabric.)") label(3 "Conflict Prob. Low Aid Prob. (fabric.)")  cols(1) size(3)) name(gr`t', replace)   graphregion( color(white) ) bgcolor(white)
graph export "$graphs\low_vs_high_brd_fabricated_China1_conflict.png", as(png) replace




***************************
* Main Results IV         *
***************************

* Figure 15: Coefficient Plot – Instrumental Variable
use "$data\working_file.dta", clear
eststo clear
cap drop std_lnaid std_lnaid_c 
egen std_lnaid=std(lnaid), mean(0)  std(1)
egen std_lnaid_c=std(lnaid_c), mean(0)  std(1)

gen l1std_lnaid=l1.std_lnaid
gen l2std_lnaid_c=l2.std_lnaid_c

global wblag1 			"l1.std_lnaid"
global chinalag2		"l2.std_lnaid_c"

matrix drop _all

// Matrices
mat C=J(7,2,0)
mat E=J(7,2,0)

local no = 1


// Create actor variables
foreach i in  t1 t2 t3g t3ng {
	gen inten1_adm1_`i'=0
	replace inten1_adm1_`i'=100 if brd_`i'_adm1>=5
}




foreach var in  inten1_adm1 D_tot_rds_adm1  D_repress_nl inten1_adm1_t1 inten1_adm1_t2 inten1_adm1_t3g inten1_adm1_t3ng {

* Model 2 - state, year and region FE
	reghdfe `var' $control2 c.l2.prob_cum ($wblag1 = c.l1.IDA_position_run#c.l2.prob_cum), absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(i.state#i.transaction_year i.ID_adm1_num) stages(first) old

	capture scalar b`no' = _b[$wblag1]					
	capture scalar s`no' = _se[$wblag1]				
	capture matrix C[`no', 1] =  b`no'				
	capture matrix C[`no', 2] =  s`no'		
	capture reghdfe `var' $control2 c.l2.prob_cum (l1.lnaid = c.l1.IDA_position_run#c.l2.prob_cum), absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(i.state#i.transaction_year i.ID_adm1_num) stages(first) old
	eststo `var'_WB2
	
	local no = `no'+1	
	}
	
		local no=1
foreach var in  inten1_adm1  inten1_adm1_t1 inten1_adm1_t2 inten1_adm1_t3g inten1_adm1_t3ng D_tot_rds_adm1  D_repress_nl  {

// China
* Model 2 - state, year and region FE

	reghdfe `var' $control2  c.l3.prob_cum_c ($chinalag2 = c.l3.factor1#c.l3.prob_cum_c), ///
	absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state#i.transaction_year i.state i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(i.state#i.transaction_year i.ID_adm1_num) stages(first) old
	
	capture scalar b`no' = _b[$chinalag2]					
	capture scalar s`no' = _se[$chinalag2]				
	capture matrix E[`no', 1] =  b`no'				
	capture matrix E[`no', 2] =  s`no'	
	capture reghdfe `var' $control2  c.l3.prob_cum_c (l2.lnaid_c = c.l3.factor1#c.l3.prob_cum_c), ///
	absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state#i.transaction_year i.state i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(i.state#i.transaction_year i.ID_adm1_num) stages(first) old
	eststo `var'_CHINA2

	local no = `no'+1
}

**************************
graph set window fontface stSans
coefplot (matrix(C[,1]), se(C[,2]) mcolor("58 85 164") msymbol(o) ciopts(lcolor("58 85 164") lwidth(thin) recast(rcap))  )	  , ///
graphregion(color(white))  xline(0, lcolor(grey) lwidth(thin) ) grid(none)  /// 
coeflabels(r1="Outright conflict (OC): >=5 BRD" ///
	r2="OC: State vs. N-state" ///
	r3="OC: N-state vs. N-state" /// 
	r4="OC: State vs. civilians"  /// 
	r5="OC: N-state vs. civilians" ///
	r6="Riots, demonstrations and strikes" ///
	r7="Non-lethal government repression" ) ///
	xlabel( , labsize(vsmall))   levels(90)  fxsize(152)  ///
headings(r1 = "{bf:Panel A: Types of Conflict}"  ///
r4 = "{bf:Panel B: Types of Actors}"  ) ///
xlabel( , labsize(vsmall))     ///
title("WB", size(small)) /// legend(on order(2 "Country-year and region fixed effects") size(vsmall) col(1) ) ///
xscale(range(-10(2)6)) xlabel(-10(2)6) 
graph save "$graphs\wb_iv.gph", replace

//			    legend(label(2 "State and year fixed effects")   label(4 "State, year and region fixed effects") size(tiny) col(1) )



coefplot (matrix(E[,1]), se(E[,2]) mcolor("212 82 154") msymbol(o) ciopts(lcolor("212 82 154") lwidth(thin) recast(rcap)) ) , ///
graphregion(color(white))  xline(0, lcolor(grey) lwidth(thin) ) grid(none)  ///
xlabel( , labsize(vsmall))  nolabels ylabel(,nolabel)  levels(90)  fxsize(80)  /// 
headings(r1 = " "  ///
r4 = " "  ) ///
title("China", size(small)) /// legend(on order(2 "Year and region FE" 4 "Country-year and region FE") size(vsmall) col(1) ) ///
xscale(range(-10(2)6)) xlabel(-10(2)6) 

graph save "$graphs\china_iv.gph", replace
graph combine "$graphs\wb_iv.gph" "$graphs\china_iv.gph"  , graphregion(color(white))
graph export "$graphs/coefplot_iv_final.pdf", as(pdf) replace 
graph export "$graphs/coefplot_iv_final.eps", replace 	
// note: the quality of the png is very bad when the code is executed. I opened the graph editor and save the graph as pdf, which has much better quality and a small file size 


* Table 13: Full IV – Conflict Actors and Types ADM1-Level

esttab  inten1_adm1_WB2  inten1_adm1_t1_WB2  inten1_adm1_t2_WB2  inten1_adm1_t3g_WB2  inten1_adm1_t3ng_WB2 D_tot_rds_adm1_WB2  D_repress_nl_WB2 using "$tables\IV_multi_Pss_wb1_diff_outcomes", replace booktabs fragment keep(L.lnaid) coeflabels(L.lnaid "$ln(World\,Bank\,Aid\,\textsubscript{t-1})$") ///
se star(* 0.10 ** 0.05 *** 0.01) b(%9.4f) compress nomtitles parentheses prehead("IV Second Stage: WB \\") stats(N idp widstat ,fmt(0 3 3)label("N" "Kleibergen-Paap underidentification test p-value" "Kleibergen-Paap weak identification F-statistic")) 	nonumbers noline nogap
				
esttab  inten1_adm1_CHINA2  inten1_adm1_t1_CHINA2 inten1_adm1_t2_CHINA2  inten1_adm1_t3g_CHINA2 inten1_adm1_t3ng_CHINA2  D_tot_rds_adm1_CHINA2  D_repress_nl_CHINA2 using "$tables\IV_multi_Pss_china1_diff_outcomes", replace booktabs fragment keep(L2.lnaid_c) coeflabels(L2.lnaid_c "$ln(Chinese\,Aid\,\textsubscript{t-2})$") ///
se star(* 0.10 ** 0.05 *** 0.01) b(%9.4f) compress nomtitles parentheses prehead("IV Second Stage: China \\") stats(N idp widstat ,fmt(0 3 3)label("N" "Kleibergen-Paap underidentification test p-value" "Kleibergen-Paap weak identification F-statistic")) nonumbers noline nogap 


***************************
* Robustness Checks IV:  *
***************************
* Leave One Out Analysis:  Figure 18/19: Robustness of first stage for World Bank / Chinese Aid - Leaving one country out
* IDA
use "$data\working_file.dta",clear
levelsof ISO3, local(countries)
foreach country of local countries {
	quietly eststo `country':	reghdfe l1.lnaid  $control c.l.prob_cum c.l1.IDA_position_run#c.l.prob_cum if ISO3!="`country'",  absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(fcy fr)
}

coefplot AGO || BDI || BEN || BFA || BWA || CAF || CIV || CMR || CMR || COD || COG || DZA, keep(cL.IDA_position_run#cL.prob_cum) yline(0) vertical bycoefs byopts(yrescale) graphregion( color(white) )  yla(-20 80) bgcolor(white)
graph export "$graphs\fsl1o\WB_cy_1.png", as(png) replace
coefplot  EGY || ERI || ETH || GAB || GHA || GIN || GMB || GNB || KEN || LBR || LBY || LSO, keep(cL.IDA_position_run#cL.prob_cum) yline(0) vertical bycoefs byopts(yrescale) graphregion( color(white) ) yla(-20 80) bgcolor(white)
graph export "$graphs\fsl1o\\WB_cy_2.png", as(png) replace
coefplot MAR || MDG || MLI || MOZ || MRT || MUS || MWI || NAM || NER || NGA || RWA || SEN, keep(cL.IDA_position_run#cL.prob_cum) yline(0) vertical bycoefs byopts(yrescale) graphregion( color(white) ) yla(-20 80)  bgcolor(white)
graph export "$graphs\fsl1o\\WB_cy_3.png", as(png) replace
coefplot SLE || SOM || TCD || TGO || TUN || TZA || UGA || ZAF || ZMB || ZWE , keep(cL.IDA_position_run#cL.prob_cum) yline(0) vertical bycoefs byopts(yrescale) graphregion( color(white) ) yla(-20 80) bgcolor(white)
graph export "$graphs\fsl1o\\WB_cy_4.png", as(png) replace


 * Bluhm et al Factor1
 levelsof ISO3, local(countries) 

 foreach country of local countries {
 	quietly eststo `country':	reghdfe l.lnaid_c $control  c.l3.prob_cum_c  c.l3.factor1#c.l3.prob_cum_c if ISO3!="`country'", absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq)  cluster(fcy fr)
 }
 coefplot AGO || BDI || BEN || BFA || BWA || CAF || CIV || CMR || CMR || COD || COG || DZA, keep(cL3.factor1#cL3.prob_cum_c) yline(0) vertical bycoefs byopts(yrescale) yla(-120 20) graphregion( color(white) ) bgcolor(white)
 graph export "$graphs\fsl1o\\CH1_cy_1.png", as(png) replace
 coefplot  EGY || ERI || ETH || GAB || GHA || GIN || GMB || GNB || KEN || LBR || LBY || LSO, keep(cL3.factor1#cL3.prob_cum_c) yline(0) vertical bycoefs byopts(yrescale) yla(-120 20) graphregion( color(white) ) bgcolor(white)
 graph export "$graphs\fsl1o\\CH1_cy_2.png", as(png) replace
 coefplot MAR || MDG || MLI || MOZ || MRT || MUS || MWI || NAM || NER || NGA || RWA || SEN, keep(cL3.factor1#cL3.prob_cum_c) yline(0) vertical bycoefs byopts(yrescale) yla(-120 20) graphregion( color(white) ) bgcolor(white)
 graph export "$graphs\fsl1o\\CH1_cy_3.png", as(png) replace
 coefplot SLE || SOM || TCD || TGO || TUN || TZA || UGA || ZAF || ZMB || ZWE , keep(cL3.factor1#cL3.prob_cum_c) yline(0) vertical bycoefs byopts(yrescale) yla(-120 20) graphregion( color(white) ) bgcolor(white)
 graph export "$graphs\fsl1o\\CH1_cy_4.png", as(png) replace


 ************************
 * Tables                *
 *************************
 * IV Main Results
 ********************
* IV Regressions   *
********************
* Bookmark 16.01.2022: Check if this is still in paper
* Define treatments in standard deviations
cap drop std_lnaid std_lnaid_c
egen std_lnaid=std(lnaid), mean(0)  std(1)
egen std_lnaid_c=std(lnaid_c), mean(0)  std(1)

global control2 droughtend_spi droughtstart_spi temp prec_gpcc c.ttime_mean#i.transaction_year  c.landarea#i.transaction_year  c.dist_capital#i.transaction_year c.mean_temp#i.transaction_year c.elevation_std#i.transaction_year c.Dborder#i.transaction_year c.Driver#i.transaction_year c.Dlake#i.transaction_year c.Docean#i.transaction_year c.lat#i.transaction_year c.lon#i.transaction_year

eststo clear
* IV World Bank IDA
			*** Cluster country-transaction_year and regional  ***
			* Column (8)
			reghdfe inten1_adm1 $control2 c.l2.prob_cum (l1.std_lnaid = c.l1.IDA_position_run#c.l2.prob_cum) ,  absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(fcy fr) stages(first) old
			eststo WB_ss
			estimates restore reghdfe_first1
			eststo WB_fs

				esttab WB_ss using "$tables\IV_WBss_conflict", replace booktabs fragment keep(L.std_lnaid) coeflabels(L.std_lnaid "$ln(World\,Bank\,Aid\,\textsubscript{t-1})$") ///
				se star(* 0.10 ** 0.05 *** 0.01) b(%9.4f) compress nomtitles parentheses prehead("IV Second Stage: World Bank \\") stats(N idp widstat ,fmt(0 3 3)label("N" "Kleibergen-Paap underidentification test p-value" "Kleibergen-Paap weak identification F-statistic")) nonumbers noline nogap

				esttab WB_fs using "$tables\IV_WBfs_conflict", replace booktabs fragment keep(cL.IDA_position_run#cL2.prob_cum) coeflabels(cL.IDA_position_run#cL2.prob_cum "$\,IDA\,Position\,\textsubscript{t-1}\,$\times$\,Cum.\,Prob\,\textsubscript{t-2}$") ///
				se star(* 0.10 ** 0.05 *** 0.01) b(%9.4f) compress nomtitles parentheses prehead("IV First stage: World Bank \\") ///
				nonumbers noline nogap 
			

			* IV Factor1 Bluhmetal2019
local FE  i.state `"i.state i.transaction_year"' `"i.state i.transaction_year i.ID_adm1_num"' `"i.state i.transaction_year i.ID_adm1_num i.continent_num#i.transaction_year"' 
local TimeTrend i.state#c.trend `" i.state#trend i.state#c.trend_sq"' `"i.state#c.trend i.state#c.trend_sq i.state#c.trend_cub"' `"i.state#c.trend i.state#c.trend_sq i.state#c.trend_cub i.ID_adm1_num#c.trend"' `"i.state#c.trend i.state#c.trend_sq i.state#c.trend_cub i.ID_adm1_num#c.trend i.ID_adm1_num#c.trend_sq"'

eststo clear
xtset ID_adm1_num transaction_year



			* Column (8)
			quietly reghdfe inten1_adm1 $control2  c.l3.prob_cum_c (L2.std_lnaid_c = c.l3.factor1#c.l3.prob_cum_c), ///
			absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq)  ///
			cluster(fcy fr) stages(first) old
			eststo CH_ss
			estimates restore reghdfe_first1
			eststo CH_fs
			
						esttab CH_ss using "$tables\IV_CHss_conflict", replace booktabs fragment keep(L2.std_lnaid_c) coeflabels(L2.std_lnaid_c "$ln(Chinese\,Aid\,\textsubscript{t-2})$") ///
				se star(* 0.10 ** 0.05 *** 0.01) b(%9.4f) compress nomtitles parentheses prehead("IV Second Stage: China \\") stats(N idp widstat ,fmt(0 3 3)label("N" "Kleibergen-Paap underidentification test p-value" "Kleibergen-Paap weak identification F-statistic")) ///
				nonumbers noline nogap 

				
				esttab CH_fs using "$tables\IV_CHfs_conflict", replace booktabs fragment keep(cL3.factor1#cL3.prob_cum_c) ///
				coeflabels(cL3.factor1#cL3.prob_cum_c "$\,Chinese\,Commodity\,\textsubscript{t-3}\,$\times$\,Cum.\,Prob\,\textsubscript{t-3}$") ///
				se star(* 0.10 ** 0.05 *** 0.01) b(%9.4f) compress nomtitles parentheses prehead("IV First stage: China \\") ///
				nonumbers noline nogap
				




***********************************************
* With third year as initial probability year  - Table 16: IV results - Initial Probability
use "$data\working_file.dta",clear
eststo clear
***********************************************
* Generate new stable probability from third year of panel
*** Cluster country-transaction_year and regional  ***

* Column (8)
ivreghdfe inten1_adm1 $control2 c.l2.prob_cum_98 (l1.std_lnaid = c.l1.IDA_position_run#c.l2.prob_cum_98) if transaction_year>=1998 , absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq)  cluster(fcy fr)
eststo WB_ini_inten1_cy3_ss
capture reghdfe l1.std_lnaid $control2 c.l2.prob_cum_98 c.l1.IDA_position_run#c.l2.prob_cum_98 if transaction_year>=1998 & inten1_adm1!=., absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq)  cluster(fcy fr)
eststo WB_ini_inten1_cy3_fs

esttab WB_ini_inten1_cy*_fs using "$tables\IV_multi_Pfs_wb_adm1_initial_prob", replace booktabs fragment keep(cL.IDA_position_run#cL2.prob_cum_98) ///
coeflabels(cL.IDA_position_run#cL2.prob_cum_98 "$\,IDA\,Position\,\textsubscript{t-1}\,$\times$\,P\,\textsubscript{98}$" ) ///
se star(* 0.10 ** 0.05 *** 0.01) b(%9.4f) compress   nomtitles parentheses prehead("IV First stage: World Bank \\") ///
nonumbers noline nogap 
esttab WB_ini_inten1_cy*_ss using "$tables\IV_multi_Pss_wb_inten1_adm1_initial_prob", replace booktabs fragment keep(L.std_lnaid) ///
coeflabels(L.std_lnaid "$ln(World\,Bank\,Aid\,\textsubscript{t-1})$")  ///
se star(* 0.10 ** 0.05 *** 0.01) b(%9.4f) compress   nomtitles parentheses prehead("IV Second Stage: World Bank \\") stats(idp widstat ,fmt(3 3)label("Kleibergen-Paap underidentification test p-value" "Kleibergen-Paap weak identification F-statistic"))  nonumbers noline nogap

* IV Chinese Commodity
* Column (8)
ivreghdfe inten1_adm1 $control2 c.l3.prob_cum_c_03 (l2.std_lnaid_c = c.l3.factor1#c.l3.prob_cum_c_03) if transaction_year>2003, absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq)  cluster(fcy fr)
eststo CHINA1_ini_inten1_cy3_ss
reghdfe l2.std_lnaid_c $control2 c.l3.prob_cum_c_03 c.l3.factor1#c.l3.prob_cum_c_03 if transaction_year>2003 & inten1_adm1!=., absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq)  cluster(fcy fr)
eststo CHINA1_ini_inten1_cy3_fs

esttab CHINA1_ini_inten1_cy*_fs using "$tables\IV_multi_Pfs_china1_adm1_initial_prob", replace booktabs fragment keep(cL3.factor1#cL3.prob_cum_c_03) ///
coeflabels(cL3.factor1#cL3.prob_cum_c_03 "$\,Chinese\,Commodity\,\textsubscript{t-3}\,$\times$\,P\,\textsubscript{03}$") ///
se star(* 0.10 ** 0.05 *** 0.01) b(%9.4f) compress   nomtitles parentheses prehead("IV First stage: China \\") ///
nonumbers noline nogap
esttab CHINA1_ini_inten1_cy*_ss using "$tables\IV_multi_Pss_china1_inten1_adm1_initial_prob", replace booktabs fragment keep(L2.std_lnaid_c) ///
coeflabels(L2.std_lnaid_c "$ln(Chinese\,Aid\,\textsubscript{t-2})$")  ///
se star(* 0.10 ** 0.05 *** 0.01) b(%9.4f) compress   nomtitles parentheses prehead("IV Second Stage: China \\") stats(idp widstat ,fmt(3 3)label("Kleibergen-Paap underidentification test p-value" "Kleibergen-Paap weak identification F-statistic")) ///
nonumbers noline nogap 



* Leverage plots -  Figure 18: High leverage regions
use "$data\working_file.dta",clear

reg l1.lnaid  $control2 c.l2.prob_cum c.l1.IDA_position_run#c.l2.prob_cum i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#c.trend i.state#c.trend_sq
predict leverage_c, leverage
predict r_c, rstudent 
*graph export "$graphs\trends\lvr_c.png", as(png) replace
reg l1.lnaid  $control2 c.l2.prob_cum c.l1.IDA_position_run#c.l2.prob_cum i.ID_adm1_num#c.trend i.ID_adm1_num i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq
predict leverage_cy, leverage
predict r_cy, rstudent
lvr2plot
graph export "$graphs\trends\lvr_cy.png", as(png) replace

gsort -leverage_cy
capture drop n1
generate n1 = _n

sort ID_adm1_num transaction_year
reg l2.lnaid_c  $control2  c.l3.prob_cum_c c.l3.factor1#c.l3.prob_cum_c i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#c.trend i.state#c.trend_sq
lvr2plot, mlabel(ADM1) msize(tiny)
lvr2plot, mlabel(transaction_year) msize(tiny)
graph export "$graphs\trends\lvr_c_c.png", as(png) replace
predict leverage_c_c, leverage
predict r_c_c, rstudent
reg l2.lnaid_c  $control2  c.l3.prob_cum_c c.l3.factor1#c.l3.prob_cum_c i.ID_adm1_num#c.trend i.ID_adm1_num i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq
lvr2plot
graph export "$graphs\trends\lvr_c_cy.png", as(png) replace
predict leverage_c_cy, leverage
predict r_c_cy, rstudent

gsort -leverage_c_cy
capture drop n2
generate n2 = _n



* Leverage plots comparable to GPSS
gsort -leverage_cy
capture drop n1
generate n1 = _n
*Auxiliary vars
tostring(transaction_year), gen(string_year)
gen ADM0_year=ISO3+"-"+ADM0+" "+string_year
gen ADM1_year=ISO3+"-"+ADM1+" "+string_year
generate pos = 3

replace pos = 1 if (ADM1_year ==ISO3+"-"+"Iburengerazuba"+" "+"2013")
replace pos = 4 if (ADM1_year ==ISO3+"-"+"Kara"+" "+"1997")
replace pos = 1 if (ADM1_year ==ISO3+"-"+"Rwanda"+" "+"1997")
replace pos = 7 if (ADM1_year ==ISO3+"-"+"Southern"+" "+"2013")
replace pos = 9 if (ADM1_year ==ISO3+"-"+"Southern"+" "+"1997")
replace pos = 1 if (ADM1_year ==ISO3+"-"+"Southern"+" "+"1998")
replace pos = 5 if (ADM1_year ==ISO3+"-"+"Southern"+" "+"2012")
replace pos = 5 if (ADM1_year ==ISO3+"-"+"Eastern"+" "+"2012")
replace pos = 12 if (ADM1_year ==ISO3+"-"+"Eastern"+" "+"1997")
replace pos = 6 if (ADM1_year ==ISO3+"-"+"Eastern"+" "+"2013")
replace pos = 2 if (ADM1_year ==ISO3+"-"+"Eastern"+" "+"1998")
replace pos = 10 if (ADM1_year ==ISO3+"-"+"Northern"+" "+"2013")
replace pos = 10 if (ADM1_year ==ISO3+"-"+"Northern"+" "+"1997")
replace pos = 11 if (ADM1_year ==ISO3+"-"+"Western"+" "+"1997")
replace pos = 6 if (ADM1_year ==ISO3+"-"+"Umujyi wa Kigali"+" "+"1997")
replace pos = 4 if (ADM1_year ==ISO3+"-"+"Umujyi wa Kigali"+" "+"2013")
replace pos = 11 if (ADM1_year ==ISO3+"-"+"Western"+" "+"2013")
replace pos = 6 if (ADM1_year ==ISO3+"-"+"Agadez"+" "+"2013")
twoway scatter r_cy leverage_cy  if n1 <= 30, jitter(10) mfcolor(blue) mlcolor(dblue) mlabel(ADM1_year) mlabgap(0.2cm) msize(0.6) mlabsize(1.7) mlabv(pos)  xtitle("Studentized Residuals", height(7) size(small)) ytitle("Leverage", height(7) size(small)) xlabel(,labsize(small)) ylabel(,labsize(small)) yline(10, lcolor(black) lpattern(dash)) legend(off) graphregion( color(white))
graph export "$graphs\leverage_top30_wb_cy_ADM1.png", as(png) replace

* China
gsort -leverage_c_cy 
capture drop n2
generate n2 = _n


replace pos = 7 if (ADM1_year ==ISO3+"-"+"Southern"+" "+"2013")
replace pos = 2 if (ADM1_year ==ISO3+"-"+"Iburengerazuba"+" "+"2013")
replace pos = 1 if (ADM1_year ==ISO3+"-"+"Plateaux"+" "+"2013")
replace pos = 4 if (ADM1_year ==ISO3+"-"+"Kara"+" "+"2013")
replace pos = 6 if (ADM1_year ==ISO3+"-"+"Kara"+" "+"2003")
replace pos = 6 if (ADM1_year ==ISO3+"-"+"Western"+" "+"2003")
replace pos = 2 if (ADM1_year ==ISO3+"-"+"Plateaux"+" "+"2003")
replace pos = 1 if (ADM1_year ==ISO3+"-"+"Rwanda"+" "++" "+"1997")
replace pos = 6 if (ADM1_year ==ISO3+"-"+"Southern"+" "+"2003")
replace pos = 4 if (ADM1_year ==ISO3+"-"+"Amajyaruguru"+" "+"2003")
replace pos = 1 if (ADM1_year ==ISO3+"-"+"Southern"+" "+"1998")
replace pos = 5 if (ADM1_year ==ISO3+"-"+"Southern"+" "+"2012")
replace pos = 12 if (ADM1_year ==ISO3+"-"+"Eastern"+" "+"2003")
replace pos = 6 if (ADM1_year ==ISO3+"-"+"Eastern"+" "+"2013")
replace pos = 5 if (ADM1_year ==ISO3+"-"+"Northern"+" "+"2003")
replace pos = 11 if (ADM1_year ==ISO3+"-"+"Northern"+" "+"2013")
replace pos = 3 if (ADM1_year ==ISO3+"-"+"Savanes"+" "+"2003")
replace pos = 5 if (ADM1_year ==ISO3+"-"+"Centre"+" "+"2003")
replace pos = 3 if (ADM1_year ==ISO3+"-"+"Umujyi wa Kigali"+" "+"2003")
replace pos = 11 if (ADM1_year ==ISO3+"-"+"Western"+" "+"2013")
replace pos = 6 if (ADM1_year ==ISO3+"-"+"Anseba"+" "+"2003")
twoway scatter r_c_cy leverage_c_cy  if n2 <= 30, mfcolor(blue) mlcolor(dblue) mlabel(ADM1_year) mlabgap(0.2cm) msize(0.6) mlabsize(1.7) mlabv(pos) xtitle("Studentized Residuals", height(7) size(small)) ytitle("Leverage", height(7) size(small)) xlabel(,labsize(small)) ylabel(,labsize(small)) yline(10, lcolor(black) lpattern(dash)) legend(off) graphregion( color(white)) 
graph export "$graphs\leverage_top30_china_cy_ADM1.png", as(png) replace
graph close
*Drop auxiliary vars
drop pos
drop string_year
drop ADM0_year
drop ADM1_year
*drop clock
*drop clock_auto



*******************************
* Without high leverage (wohl) / 
*******************************
eststo clear
xtset ID_adm1_num transaction_year
sort ID_adm1_num transaction_year
*Setting it up to exclude top 30, 1%, 2%, 5% leverage from sample
foreach i in 30 123 616 {
	* Column (6)
	/*reghdfe inten1_adm1 $control2  c.l3.prob_cum_c (l2.lnaid_c = c.l3.factor1#c.l3.prob_cum_c) if n2>`i', absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state i.transaction_year i.state#c.trend i.state#c.trend_sq) stages(first) cluster(i.state#i.transaction_year i.ID_adm1_num) old
	eststo CHINA1_`var'_cy2_ss`i'
	estimates restore reghdfe_first1
	eststo CHINA1_`var'_cy2_fs`i'*/
	* Column (8)
	reghdfe inten1_adm1 $control2 c.l2.prob_cum (l1.std_lnaid = c.l1.IDA_position_run#c.l2.prob_cum) if n1>`i', absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) stages(first) cluster(i.state#i.transaction_year i.ID_adm1_num) old
	eststo IDAP_inten1_cy3_ss`i'
	estimates restore reghdfe_first1
	eststo IDAP_inten1_cy3_fs`i'
}
foreach i in 30 79 398 {
	reghdfe inten1_adm1 $control2  c.l3.prob_cum_c (l2.std_lnaid_c = c.l3.factor1#c.l3.prob_cum_c) if n2>`i', absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) stages(first) cluster(i.state#i.transaction_year i.ID_adm1_num) old
	eststo CHINA1_`var'_cy3_ss`i'
	estimates restore reghdfe_first1
	eststo CHINA1_`var'_cy3_fs`i'
	*esttab CHINA1_`var'_cy*_fs using "$tables\IV_multi_Pfs_china1_adm1_wohl`i'", replace booktabs fragment keep(cL3.factor1#cL3.prob_cum_c L3.prob_cum_c) ///
	*coeflabels(cL3.factor1#cL3.prob_cum_c "$\,Chinese\,Commodity\,\textsubscript{t-3}\,$\times$\,Cum.\,Prob\,\textsubscript{t-3}$" L3.prob_cum_c "$\,Cum.\,Prob\,\textsubscript{t-3}$") ///
	*se star(* 0.10 ** 0.05 *** 0.01) b(%9.4f) compress   nomtitles parentheses prehead("IV First stage: China \\") ///
	*nonumbers noline nogap
}
* Table 14: ADM1 IV World Bank - Without High Leverage Regions
esttab IDAP_inten1_cy*_fs* using "$tables\IV_multi_Pfs_wb_adm1_noprob_wohl", replace booktabs fragment keep(cL.IDA_position_run#cL2.prob_cum) ///
coeflabels(cL.IDA_position_run#cL2.prob_cum "$\,IDA\,Pos\,\textsubscript{t-1}\,$\times$\,P\,\textsubscript{t-2}$") ///
se star(* 0.10 ** 0.05 *** 0.01) b(%9.4f) compress   nomtitles parentheses prehead("IV First stage\\") ///
nonumbers noline nogap 
esttab IDAP_inten1_cy*_ss* using "$tables\IV_multi_Pss_wb_inten1_adm1_wohl", replace booktabs fragment keep(L.std_lnaid) ///
coeflabels(L.std_lnaid "$ln(World\,Bank\,Aid\,\textsubscript{t-1})$") ///
se star(* 0.10 ** 0.05 *** 0.01) b(%9.4f) compress   nomtitles parentheses prehead("IV Second Stage\\") stats(idp widstat ,fmt(3 3)label("KP p-value" "KP F-stat"))  nonumbers noline nogap

* Table 15: ADM1 IV China - Without High Leverage Regions
esttab CHINA1_`var'_cy*_fs* using "$tables\IV_multi_Pfs_china1_adm1_noprob_wohl", replace booktabs fragment keep(cL3.factor1#cL3.prob_cum_c) ///
coeflabels(cL3.factor1#cL3.prob_cum_c "$\,Comm\,\textsubscript{t-3}\,$\times$\,P\,\textsubscript{t-3}$") ///
se star(* 0.10 ** 0.05 *** 0.01) b(%9.4f) compress   nomtitles parentheses prehead("IV First stage\\") ///
nonumbers noline nogap
esttab CHINA1_`var'_cy*_ss* using "$tables\IV_multi_Pss_china1_inten1_adm1_wohl", replace booktabs fragment keep(L2.std_lnaid_c) ///
coeflabels(L2.std_lnaid_c "$ln(Chinese\,Aid\,\textsubscript{t-2})$")  ///
se star(* 0.10 ** 0.05 *** 0.01) b(%9.4f) compress   nomtitles parentheses prehead("IV Second Stage\\") stats(idp widstat ,fmt(3 3)label("KP p-value" "KP F-stat")) ///
nonumbers noline nogap 

///////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////

*******************************
* Robustness Coeff Plot 	  *
*******************************
* Figure 17: Coefficient Plot First Stages – Instrumental Variable
eststo clear
*Getting prepared for future coeff plot
mata: mata clear
matrix drop _all
*WB matrixmat to collect estimates
mat W=J(10,2,0)
*CH matrix to collect estimates
mat C=J(10,2,0)
*counting local
local no = 1
*Initial probability
reghdfe l1.std_lnaid $control2 c.l2.prob_cum_98   c.l1.IDA_position_run#c.l2.prob_cum_98 if transaction_year>=1998 & inten1_adm1!=., absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq)  cluster(fcy fr)
eststo WB_ini_inten1_cy3_fs
*WB parameters of interest
capture scalar bw`no' = _b[c.l1.IDA_position_run#c.l2.prob_cum_98]
capture scalar sw`no' = _se[c.l1.IDA_position_run#c.l2.prob_cum_98]
capture matrix W[`no', 1] =  bw`no'
capture matrix W[`no', 2] =  sw`no'
reghdfe l2.std_lnaid_c $control2 c.l3.prob_cum_c_03  c.l3.factor1#c.l3.prob_cum_c_03 if transaction_year>2003 & inten1_adm1!=., absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq)  cluster(fcy fr)
eststo CHINA1_ini_inten1_cy3_fs
*CH parameters of interest
capture scalar bc`no' = _b[c.l3.factor1#c.l3.prob_cum_c_03]
capture scalar sc`no' = _se[c.l3.factor1#c.l3.prob_cum_c_03]
capture matrix C[`no', 1] =  bc`no'
capture matrix C[`no', 2] =  sc`no'
local no=`no'+1
*Leverages
foreach i in 30 123 616 {
	reghdfe inten1_adm1 $control2 c.l2.prob_cum (l1.std_lnaid = c.l1.IDA_position_run#c.l2.prob_cum) if n1>`i', absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) stages(first) cluster(i.state#i.transaction_year i.ID_adm1_num) old
	estimates restore reghdfe_first1
	eststo WB1_`var'_cy3_fs`i'
	*WB parameters of interest
	capture scalar bw`no' = _b[c.l1.IDA_position_run#c.l2.prob_cum]
	capture scalar sw`no' = _se[c.l1.IDA_position_run#c.l2.prob_cum]
	capture matrix W[`no', 1] =  bw`no'
    capture matrix W[`no', 2] =  sw`no'
	local no=`no'+1
}
local no=`no'-3
foreach i in 30 79 398 {
	reghdfe inten1_adm1 $control2  c.l3.prob_cum_c (l2.std_lnaid_c = c.l3.factor1#c.l3.prob_cum_c) if n2>`i', absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) stages(first) cluster(i.state#i.transaction_year i.ID_adm1_num) old
	estimates restore reghdfe_first1
	eststo CHINA1_`var'_cy3_fs`i'
	*CH parameters of interest
	capture scalar bc`no' = _b[c.l3.factor1#c.l3.prob_cum_c]
	capture scalar sc`no' = _se[c.l3.factor1#c.l3.prob_cum_c]
	capture matrix C[`no', 1] =  bc`no'
	capture matrix C[`no', 2] =  sc`no'
	local no=`no'+1
}
*Further Controls
capture rename global_temperature_anonmaly globaltemp
eststo clear
qui{
	foreach int in globaltemp asinhbrd_global asinhfdi_glob {
		reghdfe l1.std_lnaid $control2 c.l2.prob_cum c.l1.IDA_position_run#c.l2.prob_cum c.l2.prob_cum#c.l1.`int' ,  absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(fcy fr)
		eststo WB_fs_`int'
		capture scalar bw`no' = _b[c.l1.IDA_position_run#c.l2.prob_cum]
		capture scalar sw`no' = _se[c.l1.IDA_position_run#c.l2.prob_cum]
		capture matrix W[`no', 1] =  bw`no'
		capture matrix W[`no', 2] =  sw`no'
		reghdfe l2.std_lnaid_c $control2 c.l3.prob_cum_c   c.l3.factor1#c.l3.prob_cum_c c.l3.prob_cum_c#c.l3.`int', absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(fcy fr)
		eststo CH_fs_`int'
		capture scalar bc`no' = _b[c.l3.factor1#c.l3.prob_cum_c]
		capture scalar sc`no' = _se[c.l3.factor1#c.l3.prob_cum_c]
		capture matrix C[`no', 1] =  bc`no'
		capture matrix C[`no', 2] =  sc`no'
		local no=`no'+1
	}
	foreach i in nat_res_s lights_sum isum_pop_ADM1 {
		reghdfe l1.std_lnaid $control2 c.l2.prob_cum c.l1.IDA_position_run#c.l2.prob_cum c.l2.`i' c.l2.`i'#c.l1.IDA_position_run,  absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(fcy fr)
		eststo WB1_fs_`i'
		capture scalar bw`no' = _b[c.l1.IDA_position_run#c.l2.prob_cum]
		capture scalar sw`no' = _se[c.l1.IDA_position_run#c.l2.prob_cum]
		capture matrix W[`no', 1] =  bw`no'
		capture matrix W[`no', 2] =  sw`no'
		reghdfe l2.std_lnaid_c $control2 c.l3.prob_cum_c c.l3.factor1#c.l3.prob_cum_c c.l3.`i' c.l3.`i'#c.l3.factor1, absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(fcy fr)
		eststo CH1_fs_`i'
		capture scalar bc`no' = _b[c.l3.factor1#c.l3.prob_cum_c]
		capture scalar sc`no' = _se[c.l3.factor1#c.l3.prob_cum_c]
		capture matrix C[`no', 1] =  bc`no'
		capture matrix C[`no', 2] =  sc`no'
		local no=`no'+1
	}
}

graph set window fontface stSans

coefplot (matrix(W[,1]), se(W[,2]) mcolor("58 85 164") msymbol(o) msize(small) ciopts(lcolor("58 85 164") lwidth(thin) recast(rcap))  ) 	  , ///
graphregion(color(white))  xline(0, lcolor(grey) lwidth(vthin) ) grid(none) /// 
coeflabels(r1="Without Top 30 High Leverage" ///
	r2="Without Top 1% High Leverage" ///
	r3="Without Top 5% High Leverage" ///
	r4="Initial Probability" ///
	r5="Controlled for Probability x Global T°" ///
	r6="Controlled for Probability x Global Conflict" ///
	r7="Controlled for Probability x Global FDI" ///
	r8="Controlled for Natural Resources x Time Series" ///
	r9="Controlled for Lights x Time Series" ///
	r10="Controlled for Population x Time Series", labsize(small) ) ///
	xlabel( , labsize(small))   levels(90)  fxsize(157)  ///
headings(r1 = "{bf:Panel A: Influence of Outliers}" ///
r4 = "{bf:Panel B: Computation of Probability}" ///
r5 = "{bf:Panel C: Role of the Cross-sectional Terms}" ///
r8 = "{bf:Panel D: Role of the Exogenous Time Series}") ///
xlabel( , labsize(vsmall))     ///
title("WB", size(small)) xscale(range(-4(2)11)) xlabel(-4(2)12) 
graph save "$graphs\wb_coefplot_red.gph", replace

coefplot (matrix(C[,1]), se(C[,2]) mcolor("212 82 154") msymbol(o)  msize(small) ciopts(lcolor("212 82 154") lwidth(thin) recast(rcap)) ) , ///
graphregion(color(white))  xline(0, lcolor(grey) lwidth(thin) ) grid(none)  ///
xlabel( , labsize(small))  nolabels ylabel(,nolabel)  levels(90)  fxsize(74)  /// 
headings(r1 = "{bf:}" ///
r4 = "{bf:}" ///
r5 = "{bf:}" ///
r8 = "{bf:}") ///
xlabel( , labsize(vsmall))     ///
title("China", size(small)) xscale(range(-12(2)4)) xlabel(-12(2)4)
graph save "$graphs\china_coefplot_red.gph", replace

graph combine "$graphs/wb_coefplot_red.gph" "$graphs/china_coefplot_red.gph"  , graphregion(color(white))
graph export "$graphs/coefplot_red.pdf", as(pdf) replace 
graph export "$graphs/coefplot_red.eps", replace 	



***********************
/* Placebo             //
***********************

**# Additional first stage controls
*/
********************************************************************************
* Table 17: WB with Further Controls - Full Regression Results
capture rename global_temperature_anonmaly globaltemp
eststo clear
qui{
	foreach int in globaltemp asinhbrd_global asinhfdi_glob {
		reghdfe inten1_adm1 $control2  c.l2.prob_cum#c.l1.`int' c.l2.prob_cum  (l1.std_lnaid = c.l1.IDA_position_run#c.l2.prob_cum), absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) stages(first) cluster(i.ID_adm1_num i.transaction_year#i.state) old
		eststo WB_ss_`int'
		reghdfe l1.std_lnaid $control2 c.l2.prob_cum c.l1.IDA_position_run#c.l2.prob_cum c.l2.prob_cum#c.l1.`int' ,  absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(fcy fr)
		eststo WB_fs_`int'
		capture matrix W[`no', 1] =  bw`no'
		capture matrix W[`no', 2] =  sw`no'

		reghdfe inten1_adm1 $control2  c.l3.prob_cum_c#c.l3.`int' c.l3.prob_cum_c  (l2.std_lnaid_c = c.l3.factor1#c.l3.prob_cum_c), absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) stages(first) cluster(i.ID_adm1_num i.transaction_year#i.state) old
		eststo CH_ss_`int'
		reghdfe l2.std_lnaid_c $control2 c.l3.prob_cum_c   c.l3.factor1#c.l3.prob_cum_c c.l3.prob_cum_c#c.l3.`int', absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(fcy fr)
		eststo CH_fs_`int'
		capture matrix C[`no', 1] =  bc`no'
		capture matrix C[`no', 2] =  sc`no'
		local no=`no'+1
	}
	foreach i in nat_res_s lights_sum isum_pop_ADM1 {
		reghdfe inten1_adm1 $control2  c.l2.`i'#c.l1.IDA_position_run l2.prob_cum c.l2.`i' (l1.std_lnaid= c.l1.IDA_position_run#c.l2.prob_cum), absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) stages(first) cluster(i.ID_adm1_num i.transaction_year#i.state) old
		eststo WB1_ss_`i'
		reghdfe l1.std_lnaid $control2 c.l2.prob_cum c.l1.IDA_position_run#c.l2.prob_cum c.l2.`i' c.l2.`i'#c.l1.IDA_position_run,  absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(fcy fr)
		eststo WB1_fs_`i'
		capture matrix W[`no', 1] =  bw`no'
		capture matrix W[`no', 2] =  sw`no'

		reghdfe inten1_adm1 $control2  c.l3.`i'#c.l3.factor1 l3.prob_cum_c c.l3.`i' (l2.std_lnaid_c = c.l3.factor1#c.l3.prob_cum_c), absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) stages(first) cluster(i.ID_adm1_num i.transaction_year#i.state) old
		eststo CH1_ss_`i'
		reghdfe l2.std_lnaid_c $control2 c.l3.prob_cum_c c.l3.factor1#c.l3.prob_cum_c c.l3.`i' c.l3.`i'#c.l3.factor1, absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(fcy fr)
		eststo CH1_fs_`i'
		capture matrix C[`no', 1] =  bc`no'
		capture matrix C[`no', 2] =  sc`no'
		local no=`no'+1

	}
}
esttab WB_ss_* WB1_ss_* using "$tables/WBss_finalcontrol.tex", replace booktabs fragment keep(L.std_lnaid) coeflabels(L.std_lnaid "$ln(World\,Bank\,Aid\,\textsubscript{t-1})$") ///
se star(* 0.10 ** 0.05 *** 0.01) b(%9.3f) compress nomtitles parentheses stats(idp widstat ,fmt(3 3)label("KP p-value" "KP F-stat")) nonumbers noline nogap prehead("IV Second stage\\")
esttab WB_fs_* WB1_fs_* using "$tables/WBfs_finalcontrol.tex", replace booktabs fragment keep(cL.IDA_position_run#cL2.prob_cum) ///
coeflabels(cL.IDA_position_run#cL2.prob_cum "$\,IDA\,Pos\,\textsubscript{t-1}\,$\times$\,P\,\textsubscript{t-2}$") ///
se star(* 0.10 ** 0.05 *** 0.01) b(%9.4f) compress nomtitles parentheses prehead("IV First stage\\")  nonumbers noline nogap
esttab CH_ss_* CH1_ss_* using "$tables/CHss_finalcontrol.tex", replace booktabs fragment keep(L2.std_lnaid_c) coeflabels(L2.std_lnaid_c "$ln(Chinese\,Aid\,\textsubscript{t-2})$") ///
se star(* 0.10 ** 0.05 *** 0.01) b(%9.3f) compress nomtitles parentheses stats(idp widstat ,fmt(3 3)label("KP p-value" "KP F-stat")) nonumbers noline nogap prehead("IV Second stage\\")
esttab CH_fs_* CH1_fs_* using "$tables/CHfs_finalcontrol", replace booktabs fragment keep(cL3.factor1#cL3.prob_cum_c) ///
coeflabels(cL3.factor1#cL3.prob_cum_c "$\,Comm\,\textsubscript{t-3}\,$\times$\,P\,\textsubscript{t-3}$") ///
se star(* 0.10 ** 0.05 *** 0.01) b(%9.4f) compress nomtitles parentheses prehead("IV First stage\\") nonumbers noline nogap 

********************************************************************************
* Table 18: China with Further Controls - Full Regression Results
capture drop auxvar
capture rename steel_prod_ln_detrend steel_detrended
eststo clear
qui{
	foreach i in factor1 factor1_levs  steel_detrended us_steel_indx {
		reghdfe inten1_adm1 $control2 l3.prob_cum_c  (l2.std_lnaid_c = c.l3.`i'#c.l3.prob_cum_c), absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) stages(first) cluster(i.ID_adm1_num i.transaction_year#i.state) old
		eststo IV1_ss_`i'
		capture gen auxvar = c.l3.`i'#c.l3.prob_cum_c
		replace auxvar = c.l3.`i'#c.l3.prob_cum_c
		reghdfe l2.std_lnaid_c $control2  c.l3.prob_cum_c   auxvar, absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(fcy fr)
		eststo IV1_fs_`i'
		replace `i' = auxvar
	}
esttab IV1_ss_* using "$tables/IV_china1_SSmixture.tex", replace booktabs fragment keep(L2.std_lnaid_c) coeflabels(L2.std_lnaid_c "$ln(Chinese\,Aid\,\textsubscript{t-2})$") ///
se star(* 0.10 ** 0.05 *** 0.01) b(%9.3f) compress nomtitles parentheses stats(idp widstat ,fmt(3 3)label("KP p-value" "KP F-stat")) nonumbers noline nogap prehead("IV Second stage\\") 
esttab IV1_fs_* using "$tables/IV_china1_FSmixture.tex", replace booktabs fragment keep(auxvar) ///
coeflabels(auxvar "$\,Time\,Series\,\textsubscript{t-3}\,$\times$\,P\,\textsubscript{t-3}$") ///
se star(* 0.10 ** 0.05 *** 0.01) b(%9.4f) compress nomtitles parentheses prehead("IV First stage\\")  			nonumbers noline nogap
}

*/




**********************************************************************************************************
***Randomization China
	****  Randomization: within years, reorder districts
clear
use "$data/working_file.dta",clear
mat A3 = J(2000,1,.)
mat B3 = J(2000,1,.)
qui{
	forvalues i=1(1)2000 {
		cap drop random rindex newIV
		bys transaction_year (ID_adm1_num): gen random = runiform()
		bys transaction_year (random): gen rindex = _n
		bys transaction_year (ID_adm1_num): gen newIV = factor1*prob_cum_c[rindex]
		qui xtset
		reghdfe l2.std_lnaid_c $control2 c.l3.prob_cum_c c.l3.newIV, absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(fcy fr)
		mat A3[`i',1] = _b[c.l3.newIV]	
		mat B3[`i',1] = _se[c.l3.newIV]	
	}
}
mat li A3
mat li B3
svmat A3
svmat B3
save "$data/temp_randomization_1.dta", replace
clear
u "$data/temp_randomization_1.dta", clear
mean A3
gen t3 = abs(A3/B3)
sort ID_adm1_num transaction_year
reghdfe l2.std_lnaid_c $control2 c.l3.prob_cum_c c.l3.factor1#c.l3.prob_cum_c, absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(fcy fr)
scalar b2=_b[c.l3.factor1#c.l3.prob_cum_c]
*scalar bwow=_b[c.l3.factor1#c.l3.prob_cum_c]
*scalar swow=_se[c.l3.factor1#c.l3.prob_cum_c]
*gen p3 = (t3>  abs(bwow/swow)) if !missing(t3)
*scalar wowwow=abs(bwow/swow)
gen p3 = (t3>  abs(_b[c.l3.factor1#c.l3.prob_cum_c]/_se[c.l3.factor1#c.l3.prob_cum_c])) if !missing(t3)
/*The p-values are estimated as the proportion of times that the absolute value of the t-statistics in replication data exceed the absolute value of the original t-statistic.*/	
sum p3
gen pv3=r(mean)
scalar pv3=pv3
twoway kdensity A3, range(-3.5 0.7) title("") scheme(s2mono) graphregion(color(white)) ytitle("Density", size(mediumlarge)) xtitle("Point Coefficent", size(mediumlarge)) ///
xline(`=scalar(b2)',lpattern(dash)) note("p-value < 0.01", size(medium))
*xline(`=scalar(b2)',lpattern(dash)) note("p-value: `=scalar(pv3)'", size(mediumlarge))
graph export "$graphs/randomization_CH_prob.pdf", replace
****  Randomization: within districts, reorder yearclear
use "$data/working_file.dta",clear
mat A2 = J(5000,1,.)
mat B2 = J(5000,1,.)
qui{
	forvalues i=1(1)5000 {
		cap drop random rindex newIV
		bys ID_adm1_num (transaction_year): gen random = runiform()
		bys ID_adm1_num (random): gen rindex = _n
		bys ID_adm1_num (transaction_year): gen newIV = factor1[rindex]*prob_cum_c
		qui xtset
		reghdfe l2.std_lnaid_c $control2 c.l3.prob_cum_c c.l3.newIV, absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(fcy fr)
		mat A2[`i',1] = _b[c.l3.newIV]	
		mat B2[`i',1] = _se[c.l3.newIV]		
	}
}
mat li A2
mat li B2
svmat A2
svmat B2
save "$data/temp_randomization_2.dta", replace
clear
u "$data/temp_randomization_2.dta", clear
mean A2
gen t2 = abs(A2/B2)
reghdfe l2.std_lnaid_c $control2 c.l3.prob_cum_c c.l3.factor1#c.l3.prob_cum_c, absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(fcy fr)
scalar b2=_b[c.l3.factor1#c.l3.prob_cum_c]
gen p2 = (t2>  abs(_b[c.l3.factor1#c.l3.prob_cum_c]/_se[c.l3.factor1#c.l3.prob_cum_c])) if !missing(t2)
sum p2	/*The p-values are estimated as the proportion of times that the absolute value of the t-statistics in replication data exceed the absolute value of the original t-statistic.*/	
gen pv2=r(mean)
scalar pv2=pv2		
twoway kdensity A2, range(-3.5 0.7) title("") scheme(s2mono) graphregion(color(white)) ytitle("Density", size(mediumlarge)) xtitle("Point Coefficent", size(mediumlarge)) ///
xline(`=scalar(b2)',lpattern(dash)) note("p-value < 0.01", size(medium))
graph export "$graphs/randomization_CH_timeseries.pdf", replace

***Randomization WB
	****  Randomization: within years, reorder districts
clear
use "$data/working_file.dta",clear
mat A3 = J(5000,1,.)
mat B3 = J(5000,1,.)
sort ID_adm1_num transaction_year
gen l1_IDA_pos=c.l1.IDA_position_run
gen l2_prob_cum=c.l2.prob_cum
qui{
	forvalues i=1(1)5000 {
		cap drop random rindex newIV
		bys transaction_year (ID_adm1_num): gen random = runiform()
		bys transaction_year (random): gen rindex = _n
		bys transaction_year (ID_adm1_num): gen newIV = l1_IDA_pos*l2_prob_cum[rindex]
		qui xtset
		reghdfe l1.std_lnaid $control2 c.l3.prob_cum newIV, absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(fcy fr)
		mat A3[`i',1] = _b[newIV]	
		mat B3[`i',1] = _se[newIV]	
	}
}
mat li A3
mat li B3
svmat A3
svmat B3
save "$data/temp_randomization_3.dta", replace
clear
u "$data/temp_randomization_3.dta", clear
mean A3
gen t3 = abs(A3/B3)
reghdfe l1.std_lnaid $control2 c.l3.prob_cum c.l1.IDA_position_run#c.l2.prob_cum, absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(fcy fr)
scalar b2=_b[c.l1.IDA_position_run#c.l2.prob_cum]
gen p3 = (t3>  abs(_b[c.l1.IDA_position_run#c.l2.prob_cum]/_se[c.l1.IDA_position_run#c.l2.prob_cum])) if !missing(t3)
/*The p-values are estimated as the proportion of times that the absolute value of the t-statistics in replication data exceed the absolute value of the original t-statistic.*/	
sum p3
gen pv3=r(mean)
scalar pv3=pv3
twoway kdensity A3, range(-0.7 3) title("") scheme(s2mono) graphregion(color(white)) ytitle("Density", size(mediumlarge)) xtitle("Point Coefficent", size(mediumlarge)) ///
xline(`=scalar(b2)',lpattern(dash)) note("p-value < 0.01", size(medium))
graph export "$graphs/randomization_WB_prob.pdf", replace
****  Randomization: within districts, reorder years
clear
use "$data/working_file.dta",clear
mat A2 = J(5000,1,.)
mat B2 = J(5000,1,.)
sort ID_adm1_num transaction_year
gen l1_IDA_pos=c.l1.IDA_position_run
gen l2_prob_cum=c.l2.prob_cum
qui{
	forvalues i=1(1)5000 {
		cap drop random rindex newIV
		bys ID_adm1_num (transaction_year): gen random = runiform()
		bys ID_adm1_num (random): gen rindex = _n
		bys ID_adm1_num (transaction_year): gen newIV = l1_IDA_pos[rindex]*l2_prob_cum
		qui xtset
		reghdfe l1.std_lnaid $control2 c.l2.prob_cum newIV, absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(fcy fr)
		mat A2[`i',1] = _b[newIV]	
		mat B2[`i',1] = _se[newIV]		
	}
}
mat li A2
mat li B2
svmat A2
svmat B2
save "$data/temp_randomization_4.dta", replace
clear
u "$data/temp_randomization_4.dta", clear
mean A2
gen t2 = abs(A2/B2)
reghdfe l1.std_lnaid $control2 c.l3.prob_cum c.l1.IDA_position_run#c.l2.prob_cum, absorb(i.ID_adm1_num#c.trend i.ID_adm1_num i.state i.transaction_year i.state#i.transaction_year i.state#c.trend i.state#c.trend_sq) cluster(fcy fr)
scalar b2=_b[c.l1.IDA_position_run#c.l2.prob_cum]
gen p2 = (t2>  abs(_b[c.l1.IDA_position_run#c.l2.prob_cum]/_se[c.l1.IDA_position_run#c.l2.prob_cum])) if !missing(t2)
sum p2	/*The p-values are estimated as the proportion of times that the absolute value of the t-statistics in replication data exceed the absolute value of the original t-statistic.*/	
gen pv2=r(mean)
scalar pv2=pv2			
twoway kdensity A2, range(-0.7 3) title("") scheme(s2mono) graphregion(color(white)) ytitle("Density", size(mediumlarge)) xtitle("Point Coefficent", size(mediumlarge)) ///
xline(`=scalar(b2)',lpattern(dash)) note("p-value < 0.01", size(medium))
graph export "$graphs/randomization_WB_timeseries.pdf", replace
